#ifndef _QGCM_NETCDFAUX_SRC_
#define _QGCM_NETCDFAUX_SRC_
#include "qgcmio.h"
#include "netcdfaux.h"

REAL getmagorder(size_t size){
	REAL rvalue = 1.0;
	size = size/10;
	while( size != 0 ){
		rvalue *= 10;
		size=size/10;
	}
	return rvalue;
};

// ................................................................................ //
const REAL & getsetmissing_val(REAL newvalue, bool setit){
	static REAL missing_val = -999; // default value
	if(setit){
		missing_val = static_cast<REAL>(newvalue);
	}
	return missing_val;
}
// ................................................................................ //
// Function creates a variable if it does not exist, returns the existing variable
// if it exists with the right attributes, dims, etc, or returns NULL if there is
// some sort of error
NcVar * ncgetaddvar(NcFile * ncfile, const char * varname, NcType type, int numdims, const NcDim ** dims){
	unsigned int numvars = ncfile->num_vars();
	for(unsigned int i = 0; i < numvars; i++){
		NcVar * var_iter = ncfile->get_var(i);
		if( strcmp(var_iter->name(), varname) == 0 ){ // may have a match
			if(var_iter->num_dims() != numdims) return NULL; // error variable mismatch
			for(int k = 0; k < var_iter->num_dims(); k++){
				if(strcmp(var_iter->get_dim(k)->name(), dims[k]->name())!= 0) return NULL;
				if(var_iter->get_dim(k)->is_unlimited() != dims[k]->is_unlimited()) return NULL;
				if(var_iter->get_dim(k)->size() != dims[k]->size()) return NULL;
			} 
			return var_iter; 
		}
	}
	//variable not found, so add it 
	return ncfile->add_var(varname, type, numdims, dims);
}
// ................................................................................ //


// ***************************************************************************************

// ---------------------------------------------------------------------------------------
// This functions below are responsible for saving the static coefficients in netcdf format.
// and or reading the static coefficients from an existing netcdf file
// It uses an nchandler to achieve this . . . .
// ---------------------------------------------------------------------------------------
template <class T, int S>
bool writefielddata(const rank<INDEX,S> * map,
		const vector< const vector<T> * > & datalist,  const int size, NcFile * ncfile,
		const char * varname, 
		const vector< pair<string, string> > & string_attributes 
		){
	if( size == 0 ) return true;
	const int edge_size = S + datalist.size();
	size_t unwrap_size = size*edge_size;
	T * unwrapped = new T[ unwrap_size ];
	size_t iter = 0;
	size_t lsize = 0;
	for(size_t i = 0; i < static_cast<size_t>(size); i++){
		bool nonzero = false;
		for(size_t x = 0; x < datalist.size(); x++){
			if((*datalist[x])[i] != 0 ){
				nonzero = true;
				break;
			}
		}
		if( nonzero ){ // write the data
			lsize++;
			const rank<INDEX,S> & rval = map[i];
			for(int pos = 0; pos < S; pos++){
				REAL fval = static_cast<REAL>( rval[pos] );
				if( fval < -126 || fval > 127 ){
					fatal_stop("Index value exceeded char limits while writing static map. Please recode to fix problem.");
				}
				//char cval = static_cast<char>( rval[pos] );
				assert( iter < unwrap_size );
				unwrapped[iter] = fval;
				++iter;
			}
			for(size_t x = 0; x < datalist.size(); x++){
				unwrapped[iter] = (*datalist[x])[i]; iter++;
			}
			// now write the values of the static weights
		}
	} 
	if(lsize == 0){
		fprintf(stderr, "Error: Length of array '%s' is zero.\n", varname);
		return false;
	}
	NcDim * dimdatalength = ncfile->add_dim( (string(varname) + string("_length")).c_str(), lsize);
	NcDim * dimdatawidth = ncfile->add_dim( (string(varname) + string("_width")).c_str(), edge_size);
	const NcDim ** dimlist = new const NcDim * [ 2 ];
	dimlist[0] = dimdatalength;
	dimlist[1] = dimdatawidth;
	NcVar * var = ncgetaddvar(ncfile, varname, ncREAL, 2, dimlist);
	delete [] dimlist;
	assert(var!=NULL);
	for(unsigned int x = 0; x < string_attributes.size(); x++){
		var->add_att( string_attributes[x].first.c_str(), string_attributes[x].second.c_str() );
	}
	if(!var->put( unwrapped, lsize, edge_size, 0, 0, 0 )){
		fatal_stop("Bad dimensions given to function NcVar::put(...)");
	}
	delete [] unwrapped;
	return true;
}; 

bool writestaticweightstofile(const bool nonlinear, const static_container * Rsc, const static_container * Ssc, 
		const map_container<REAL, 6> * Amap,
		const map_container<REAL, 3> * Qmap, const map_container<REAL, 3> * QImap,
		NcFile * ncfile)
{
	const qgcomplex * Rs = Rsc->getcoefficients();
	const rank11 * Rindexmap = Rsc->get_index_map();
	const int Rsize = Rsc->get_size();

	const qgcomplex * Ss = Ssc->getcoefficients();
	const rank11 * Sindexmap = Ssc->get_index_map();
	const int Ssize = Ssc->get_size();

	vector<REAL> dlist1( Rsize, 0 );
	vector<REAL> dlist2( Rsize, 0 );
	vector< const vector<REAL> * > dlist;
	dlist.push_back(&dlist1);
	dlist.push_back(&dlist2);
	// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	vector< pair<string, string> > string_attributes;
	string_attributes.push_back( 
			pair<string,string>( "Description", "Nonzero elements (rank11 pos, value) of the R static coefficients." ) );
	// call function twice
	for(size_t i = 0; i < static_cast<size_t>(Rsize); i++){
		dlist1[i] = Rs[i].real();
		dlist2[i] = Rs[i].imag();
	}
	if(nonlinear){
		if(! writefielddata(Rindexmap, dlist, Rsize, ncfile, "R_static_weights", string_attributes)){
			fprintf(stderr, "Warning: Error writing Rs static weights.");
		}
	}
	// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
	dlist1.resize( Ssize, 0 );
	dlist2.resize( Ssize, 0 );
	string_attributes.clear();
	string_attributes.push_back( pair<string,string>( "Description", "Nonzero elements (rank11 pos, value) of the S static coefficients." ) );
	for(size_t i = 0; i < static_cast<size_t>(Ssize); i++){
		dlist1[i] = Ss[i].real();
		dlist2[i] = Ss[i].imag();
	}
	if(nonlinear){
		if(! writefielddata(Sindexmap, dlist, Ssize, ncfile, "S_static_weights", string_attributes)){
			fprintf(stderr, "Warning: Error writing Ss static weights.");
		}
	}
	// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
	const int Asize = Amap->get_size();
	const rank6 * Aindexmap = Amap->get_index_map();
	const REAL * Adata = Amap->getcoefficients();
	dlist1.resize( Asize, 0 );
	dlist.clear();
	dlist.push_back( &dlist1 );
	for(size_t i = 0; i < static_cast<size_t>(Asize); i++){
		dlist1[i] = Adata[i];
	} 
	string_attributes.clear();
	string_attributes.push_back( 
			pair<string,string>( "Description", "Nonzero elements (rank6 pos, value) of the A operator." ) );
	if( ! writefielddata(Aindexmap, dlist, Asize, ncfile, "A_operator", string_attributes)){
		fatal_stop("Error writing A_operator static weights.");
	}
	// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
	const int Qsize = Qmap->get_size();
	const rank3 * Qindexmap = Qmap->get_index_map();
	const REAL * Qdata = Qmap->getcoefficients();
	dlist1.resize( Qsize, 0 );
	dlist.clear();
	dlist.push_back( &dlist1 );
	for(size_t i = 0; i < static_cast<size_t>(Qsize); i++){
		dlist1[i] = Qdata[i];
	} 
	string_attributes.clear();
	string_attributes.push_back( 
			pair<string,string>( "Description", "Nonzero elements (rank3 pos, value) of the Q operator." ) );
	if( ! writefielddata(Qindexmap, dlist, Qsize, ncfile, "Q_operator", string_attributes)){
		fatal_stop("Error writing Q_operator static weights.");
	}
	// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
	const int QIsize = QImap->get_size();
	const rank3 * QIindexmap = QImap->get_index_map();
	const REAL * QIdata = QImap->getcoefficients();
	dlist1.resize( QIsize, 0 );
	dlist.clear();
	dlist.push_back( &dlist1 );
	for(size_t i = 0; i < static_cast<size_t>(QIsize); i++){
		dlist1[i] = QIdata[i];
	} 
	string_attributes.clear();
	string_attributes.push_back( 
			pair<string,string>( "Description", "Nonzero elements (rank3 pos, value) of the Q operator." ) );
	if( ! writefielddata(QIindexmap, dlist, QIsize, ncfile, "QI_operator", string_attributes)){
		fatal_stop("Error writing QI_operator static weights.");
	}
	return true;
};


template <class TYPE>
NcVar * writeprofiletostaticfile(const vector<TYPE> & data, const char * name, NcFile * ncfile, NcDim * dim, NcType nctype,
		const vector< pair<string, string> > & string_attributes,
		const vector< pair<string, int> > & int_attributes
		){
	if(data.size() == 0) return NULL;
	const NcDim ** dimlist = new const NcDim * [ 1 ];
	dimlist[0] = dim;
	NcVar * var = ncgetaddvar(ncfile, name, nctype, 1, dimlist);
	delete [] dimlist;
	if(var == NULL){
		throw cgcmexception("Fatal Error: Cannot prepare nc static coeffs output file for writing.", -55);
		return var;
	}
	var->add_att("missing_value", getsetmissing_val(0, false) );
	for(unsigned int x = 0; x < string_attributes.size(); x++){
		var->add_att( string_attributes[x].first.c_str(), string_attributes[x].second.c_str() );
	}
	for(unsigned int x = 0; x < int_attributes.size(); x++){
		var->add_att( int_attributes[x].first.c_str(), int_attributes[x].second );
	}
	TYPE * farray = new TYPE[data.size()];
	for(unsigned int i = 0; i < data.size(); i++){
		farray[i] = static_cast<TYPE>(data[i]);
	}
	var->put(farray, data.size(), 0, 0, 0, 0);
	delete [] farray;
	return var;
}
// ---------------------------------------------------------------------------------------
//template void writeprofiletostaticfile(const vector<double> & data, const char *, NcFile*, NcDim*, NcType); // force explicit instantiation of double
// ---------------------------------------------------------------------------------------
// opens the file if it exists, otherwise creates a new one
double nc_static_open(model_class & mo, nchandler & nc_container, const char * fname, bool open, bool overwrite){
	// ' ' ' ' ' ' ' ' ' '  ' ' ' ' ' get global parameters/profiles 
	const unsigned int numfields = mo.get_global_param_by_name(gdb_str_NUMFIELDS).getiValue();
	const vector<REAL> & kdomain = mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN).getData();
	const vector<REAL> & dkdomain= mo.get_global_prof_by_name(gdb_prof_str_DK_DOMAIN).getData();
	const vector<REAL> & kdomaincentered = mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN_CENTERED).getData();
	const vector<int> & nkdomain = mo.get_global_prof_by_name(gdb_prof_str_NK_DOMAIN).getiData();
	const vector<REAL> & fcoeffs  = mo.get_global_prof_by_name(gdb_prof_str_FCOEFFS).getData();
	const unsigned int ksize = kdomaincentered.size();
	const unsigned int nksize= nkdomain.size();

	double nchandle = 0;

	NcFile::FileMode fmode = NcFile::Write;
	if(open){ // we are opening a static coeffs file for reading, need to check consistency 
		fmode = NcFile::ReadOnly;
		nchandle = nc_container.createncfile(fname, fmode);
		if(nchandle == 0) return 0;
	}else{ // we are writing to one
		nchandle = 0;
		if(!overwrite) { // make sure the file doesn't exist first
			nchandle = nc_container.createncfile(fname, NcFile::Write);
			if(nchandle != 0){ // file already existed error;
				return false;
			}
			nc_container.deletencfile(static_cast<int>(nchandle));
		}
		nchandle = nc_container.createncfile(fname, NcFile::Replace);
		if(nchandle == 0) return 0;
		NcFile * ncfile = nc_container.getncfile(static_cast<int>(nchandle));
		//NcDim * freeindex =  
		ncfile->add_dim("freeindex");
		NcDim * Findex = ncfile->add_dim("Findex", 4);
		//NcDim * fieldi = 
		ncfile->add_dim("fieldi", numfields);
		//NcDim * fieldj = 
		ncfile->add_dim("fieldj", numfields);
		NcDim * Kindex = ncfile->add_dim("Kindex", ksize+1);
		NcDim * kindex = ncfile->add_dim("kindex", ksize);
		NcDim * nkindex= ncfile->add_dim("nkindex", nksize);
		// NcDim * r11index =
		vector< pair<string, string> > string_attributes; // not used here but used to define attributes for the variable
		vector< pair<string, int> > int_attributes; // not used here but used to define attributes for the variable
		writeprofiletostaticfile(kdomain, gdb_prof_str_K_DOMAIN, ncfile, Kindex,ncREAL,string_attributes,int_attributes);
		writeprofiletostaticfile(dkdomain, gdb_prof_str_DK_DOMAIN, ncfile, kindex,ncREAL,string_attributes,int_attributes);
		writeprofiletostaticfile(kdomaincentered, gdb_prof_str_K_DOMAIN_CENTERED, ncfile, kindex,ncREAL,string_attributes,int_attributes);
		writeprofiletostaticfile(nkdomain, gdb_prof_str_NK_DOMAIN, ncfile, nkindex,ncREAL,string_attributes,int_attributes);   
		writeprofiletostaticfile(fcoeffs, gdb_prof_str_FCOEFFS, ncfile, Findex,ncREAL,string_attributes,int_attributes);  
		// add the fourier transform coefficient attribute
		double ft_coeff = static_cast<double>( mo.get_global_param_by_name(gdb_str_FOURIER_TRANSFORM_COEFFICIENT).getfValue());
		ncfile->add_att(gdb_str_FOURIER_TRANSFORM_COEFFICIENT, ft_coeff );
		ncfile->sync(); 
	}
	return nchandle;
};
// ---------------------------------------------------------------------------------------
unsigned int findr3element(const vector<rank3> & list, const rank3 & pos){
	static unsigned int el = 0;
	while(true){
		if(list[el] == pos) return el;
		el++;
		if(el >= list.size()) el = 0;
	}
	return 0;
};
// ---------------------------------------------------------------------------------------
bool does_global_attribute_exist( NcFile * ncfile, const char * att_name){
	const int natts = ncfile->num_atts();
	for(int i = 0; i < natts; i++){
		NcAtt * att = ncfile->get_att(i);
		if( string( att->name() ) == string( att_name ) ) return true;
	}
	return false;
};
// ---------------------------------------------------------------------------------------
bool does_dim_exist( NcFile * ncfile, const char * dim_name){
	int ndims = ncfile->num_dims();
	for(int i = 0; i < ndims; i++){
		NcDim * dim = ncfile->get_dim(i);
		if( string( dim->name() ) == string( dim_name ) ) return true;
	}
	return false;
}
// ---------------------------------------------------------------------------------------
bool nc_write_static_coeffs
(
 model_class & mo,
 const char * fname,
 bool overwrite,
 const static_container * Rsc, const static_container * Ssc,
 const map_container<REAL, 6> * Amap,
 const map_container<REAL, 3> * Qmap,
 const map_container<REAL, 3> * QImap
 ){
	//const qgcomplex * Rs, const qgcomplex * Ss, const rank11 * RSindexmap, const unsigned int & map_size)
	bool nonlinear = mo.get_global_param_by_name(gdb_str_NONLINEAR).getiValue() != 0;
	nchandler nc_container;
	double nchandle = nc_static_open( mo, nc_container, fname,  false, overwrite);
	if(nchandle == 0){
		fprintf(stderr, "Error: The static coefficients output netcdf file '%s' already exists.\n", 
				fname);
		fprintf(stderr, "Cannot overwrite existing static coefficients output file.\n");
		fprintf(stderr, "Please rename or delete file or use 'mode=overwrite' switch.\n");
		return false;
	}
	NcFile * ncfile = nc_container.getncfile(static_cast<int>(nchandle));
	// add information for the index map .....  
	const vector<rank3> & aindexmap = mo.getAfieldindices();
	if(aindexmap.size() == 0){
		fprintf(stderr, "Warning: A operator has no nonzero sets of field indices.\n");
		fprintf(stderr, "There is nonlinear data to write out.\n");
	}else{
		NcDim * indexmap = ncfile->add_dim("Amapindex", aindexmap.size());
		// conver the index map to a float
		vector<REAL> aindexmap_fl;
		for(unsigned int i = 0; i < aindexmap.size(); i++){
			const rank3 & r = aindexmap[i];
			aindexmap_fl.push_back( static_cast<REAL>(r[0])*100.0 + static_cast<REAL>(r[1])*10.0 + static_cast<REAL>(r[2]) );
		}
		vector< pair<string, string> > string_attributes;
		vector< pair<string, int> > int_attributes;
		//    string_attributes.push_back( pair<string,string>( "Description", 
		//	"A operator field index map (i.e. i,j,k unique field index triplets for which A is nonzero)." ) );
		writeprofiletostaticfile(aindexmap_fl, "Afieldindexmap", ncfile, indexmap, ncREAL, string_attributes, int_attributes);
		aindexmap_fl.clear();
	}
	// ......................................
	if(! writestaticweightstofile( nonlinear, Rsc, Ssc, Amap, Qmap, QImap, ncfile) ){
		fprintf(stderr, "Error writing static weights to file.\n"); 
		return false;
	}
	ncfile->sync();
	// finally write the values of the static coefficients to  the file ....
	ncfile->close();
	return true;
};
// ---------------------------------------------------------------------------------------
template <class T>
bool check_prof_consistency( NcFile * ncfile, const char * pname, const char * dimname, const vector<T> & prof){
	// make sure profile exists and has right number of elements and size and elements are equal
	NcVar * var = ncfile->get_var(pname);
	if(var == NULL || ! var->is_valid()){
		fprintf(stderr, "Error: Static coeffs file missing necessary profile '%s'.\n", pname);
		return false;
	}
	if( var->num_dims() != 1 ){
		fprintf(stderr, "Error: Expected 1-dimensional array, but go '%i' dimensional array.\n", var->num_dims() );
		fprintf(stderr, "Error: Static coeffs file has badly formated profile '%s'.\n", pname); return false;
	}
	if( var->type() != ncREAL ){
		fprintf(stderr, "Error: Expected variable of type 'REAL', but got something else.\n");
		fprintf(stderr, "Error: Static coeffs file contains incorrect data type for profile '%s'.\n", pname); return false;
	}
	NcDim * dim = var->get_dim(0);
	if( string(dimname) != string(dim->name() ) ){
		fprintf(stderr, "Error: Profile '%s' expects dependency on dimensional '%s' but got '%s'.\n",
				pname, dimname, var->get_dim(0)->name() );
		return false;
	}
	const unsigned int dsize = var->num_vals();
	if(dsize != prof.size() ){
		fprintf(stderr, "Error: Mismatch between current model instantiation and saved data in profile '%s'\n",
				pname);
		return false;
	}
	REAL * vals = new REAL [ dsize ];
	var->get( vals, dsize, 0, 0, 0, 0 ) ;
	for(unsigned int i = 0; i < dsize; i ++){
		if( static_cast<T>( vals[i] ) != prof[i] ){
			fprintf(stderr, "Error: Element '%i' do not agree.\n", i+1);
			fprintf(stderr, "Error: Mismatch between current model instantiation and saved data in profile '%s'\n",
					pname); return false;
		}
	}
	delete [] vals;
	// all seems good
	return true;
};
// ---------------------------------------------------------------------------------------
bool check_dim_consistency( NcFile * ncfile, const vector< pair<string, int> > & required_values ){
	for(unsigned int i = 0; i < required_values.size(); i++){
		NcDim * dim = ncfile->get_dim( required_values[i].first.c_str() );
		if( dim == NULL || ! dim->is_valid()){
			fprintf(stderr, "Error: Static coefficients data file is missing necessary dimension: '%s'.\n", 
					required_values[i].first.c_str() );
			return false;
		}
		const int rsize = required_values[i].second;
		if(rsize == 0){
			if( ! dim->is_unlimited() ){
				fprintf(stderr, "Error: Dimension '%s' in static coefficients file must be unlimited, but instead is fixed size.\n",
						required_values[i].first.c_str() );
				return false;
			}
		}else if(rsize > 0){ // fixed size
			if( dim->size() != rsize ){
				fprintf(stderr, "Error: Dimension '%s' in static coefficients file has size '%li' but current model instantiation requires size '%i'\n",
						required_values[i].first.c_str(), dim->size(), rsize );
				return false;
			}
		} // else don't care about the size at all
	}
	return true; // no problems found;
};
// ---------------------------------------------------------------------------------------
template <class T>
bool get_static_var( NcFile * ncfile, T * & data, const int required_size, const char * name, NcDim * rdim, NcType type){
	NcVar * var = ncfile->get_var(name);
	if( var == NULL || ! ( var->is_valid() ) ){
		fprintf(stderr, "Error: Variable '%s' not found in file.\n", name);
		return false;
	}
	if( var->num_vals() != required_size ){
		fprintf(stderr, "Error: Inconsistent data for variable '%s'.\n", name);
		return false;
	}
	if( string(var->get_dim(0)->name()) != string(rdim->name()) ){
		fprintf(stderr, "Error: Dimension mismatch in variable '%s'.\n", name);
		return false;
	}
	if( var->type() != type ){
		fprintf(stderr, "Error: Mismatched type while reading variable '%s'.\n", name);
		return false;
	}
	// now read the data from the file into data
	data = new T[ required_size ];
	if( ! var->get( data, required_size, 0, 0, 0, 0 ) ){
		fprintf(stderr, "Error reading variable '%s'.\n", name);
		return false;
	}
	return true;
}
// ---------------------------------------------------------------------------------------
bool open_static_coeffs_file(model_class & mo, nchandler & nc_container, double & nchandle, const char * filename){
	nchandle = nc_static_open(mo, nc_container, filename, true, false); // open file for reading ... 
	if(nchandle == 0){
		fprintf(stderr, "Cannot open '%s' for reading because the file does not exist.\n", filename);
		return false;
	}
	NcFile * ncfile = nc_container.getncfile(static_cast<int>(nchandle));
	// first check to make sure there is agreement between the global attribute FOURIER_TRANSFORM_COEFFICIENT
	// and that defined by the model here.
	if( does_global_attribute_exist( ncfile, gdb_str_FOURIER_TRANSFORM_COEFFICIENT ) ){
		double ft_coeff = static_cast<double>( mo.get_global_param_by_name(gdb_str_FOURIER_TRANSFORM_COEFFICIENT).getfValue());
		double static_value = 0;
		NcAtt * att = ncfile->get_att(gdb_str_FOURIER_TRANSFORM_COEFFICIENT);
		bool status = true;
		while(true){
			if(att->type() != ncDouble){ status=false; break; };
			if(att->num_vals() != 1){ status = false; break; };
			static_value = att->as_double(0);
			if(static_value != ft_coeff){ status = false; break; }
			break;
		}
		if(!status){
			fprintf(stderr, "Error: Mismatch between current model instantiation and static coefficients input file '%s' in the parameter '%s'.\n", filename, gdb_str_FOURIER_TRANSFORM_COEFFICIENT );
			fprintf(stderr, "The static file contains a value of '%f' while the model's value is '%f'.\n", static_value, ft_coeff);
			fprintf(stderr, "You may be inadvertently mixing fourier transform conventions.\n");
			return false;
		}
	};
	// for this function to work correctly, there are certain profiles and parameters that must match up between the file
	// and what is currently residing in the model parameter space. Namely, we need the following profiles to match:
	// K_DOMAIN, K_DOMAIN_CENTERED, K_DOMAIN_SQUARED, NK_DOMAIN, F_COEFFS, 
	// number of fields, Kindex, kindex, nkindex, Amapindex
	// check to make sure the Aindexmap is matching ...
	const int numfields = mo.get_global_param_by_name(gdb_str_NUMFIELDS).getiValue();
	const vector<REAL> &kdomain = mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN).getData();
	const vector<REAL> &dkdomain = mo.get_global_prof_by_name(gdb_prof_str_DK_DOMAIN).getData();
	//const vector<REAL> &k2 = mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN_CENTERED_SQUARED).getData();
	const vector<REAL> &kdomain_centered = mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN_CENTERED).getData();
	const vector<int> &nkdomain = mo.get_global_prof_by_name(gdb_prof_str_NK_DOMAIN).getiData();
	const vector<REAL> &fcoeffs = mo.get_global_prof_by_name(gdb_prof_str_FCOEFFS).getData();

	vector< pair< string, int > > required_values; // name, size of dim required, 0 if free, -1 if it doesn't matter
	required_values.push_back( pair< string, int >( string("Findex"), 4 ) );
	required_values.push_back( pair< string, int >( string("fieldi"), numfields ) );
	required_values.push_back( pair< string, int >( string("fieldj"), numfields ) );
	required_values.push_back( pair< string, int >( string("Kindex"), kdomain.size() ) );
	required_values.push_back( pair< string, int >( string("kindex"), kdomain_centered.size() ) );
	required_values.push_back( pair< string, int >( string("nkindex"), nkdomain.size() ) );

	if( ! check_dim_consistency( ncfile, required_values ) ){
		fprintf(stderr, "Error: Dimension mismatch between current model instantiation and static coefficients input file '%s'\n",
				filename);
		return false;
	}
	// now make sure that the profiles have values that match the current model instantantiation ...
	bool profcheck;
	profcheck = check_prof_consistency(ncfile, gdb_prof_str_K_DOMAIN, "Kindex", kdomain);
	profcheck = check_prof_consistency(ncfile, gdb_prof_str_DK_DOMAIN, "kindex", dkdomain);
	profcheck = profcheck && check_prof_consistency(ncfile, gdb_prof_str_K_DOMAIN_CENTERED, "kindex", kdomain_centered);
	profcheck = profcheck && check_prof_consistency(ncfile, gdb_prof_str_NK_DOMAIN, "nkindex", nkdomain);
	profcheck = profcheck && check_prof_consistency(ncfile, gdb_prof_str_FCOEFFS, "Findex", fcoeffs);
	if(!profcheck){
		fprintf(stderr, "Error: Bad or inconsistent data in static coefficients input file '%s'.\n", filename); return false;
	}
	// need to read the Aindexmap and store it

	if(does_dim_exist( ncfile, "Amapindex" )){ // then load the A coefficients 
		double * map;
		NcDim * sdim = ncfile->get_dim( "Amapindex" ); // existence has already been checked for
		size_t map_size = sdim->size();
		if(
				!get_static_var( ncfile, map, map_size, "Afieldindexmap", sdim, ncREAL )
		  ){
			fprintf(stderr, "Error reading static coefficients data.\n");
			return false;
		}
		vector<rank3> Aindexmap( map_size, rank3('*','*','*'));
		for(int iter = 0; iter < static_cast<int>(map_size); iter++){
			rank3 r;
			int val = static_cast<int>( map[iter] );
			for( int x = 0; x < 3; x++){
				r[x] = val%10; 
				val = val/10; 
			} 
			Aindexmap[iter] = r;
		}
		mo.setAfieldindices( Aindexmap ); // assign the nonzero index map to the model for general purposes
	};
	return true;
};
// ---------------------------------------------------------------------------------------
bool read_static_var_dims( const int WIDTH, int & width, int & length, const int rwidth, NcFile * ncfile, const char * varname){
	NcDim * sdimlength = ncfile->get_dim( string( string(varname)+"_length").c_str() ); 
	// existence has already been checked for
	NcDim * sdimwidth = ncfile->get_dim( string( string(varname)+"_width").c_str() ); 
	// existence has already been checked for
	length = sdimlength->size();
	width = sdimwidth->size();
	if( rwidth + WIDTH != width ){ // we have a mismatch
		fprintf(stderr, "Variable '%s' expects a width of '%i' but read in a width of '%i'.", varname, rwidth + WIDTH, width);
		return false;
	}
	if( length <= 0 ){
		fprintf(stderr, "Variable '%s' is empty.\n", varname);
		return false;
	}
	return true;
};

// ---------------------------------------------------------------------------------------
bool read_static_var(const char * varname, static_container & target, NcFile * ncfile){
	const int WIDTH = 11; // rank size
	const int rwidth = 2; // require 2 elements to form a complex number
	int width, length = 0;
	if(! read_static_var_dims(WIDTH, width, length, rwidth, ncfile, varname) ) return false;
	// now read in variable
	NcVar * var = ncfile->get_var(varname);
	if( var == NULL || ! ( var->is_valid() ) ){
		fprintf(stderr, "Error: Variable '%s' not found in file.\n", varname);
		return false;
	}
	vector<rank11> map( length, rank11(0,0,0,0,0,0,0,0,0,0,0) );
	vector< qgcomplex > data( length, qgcomplex(0,0) );
	int dsize = length*width;
	REAL * fdata = new REAL[dsize];
	if( ! var->get( fdata, length, width, 0, 0, 0 ) ){
		fprintf(stderr, "Error reading data from variable '%s'.\n", varname);
		delete [] fdata;
		return false;
	}
	// read data into index map etc ..
	int iter = 0, pos = 0;
	while(pos < dsize){
		rank11 r11;
		for(int i = 0; i < WIDTH; i++){
			const REAL & cval = fdata[pos];
			if(cval < -127 || cval > 128){
				fprintf(stderr, "Error reading index map element '%i': Value exceeds INDEX limits.\n", iter);
				fprintf(stderr, "Error reading data from variable '%s'.\n", varname); delete [] fdata;
				return false;
			}
			r11[i] = static_cast<INDEX>( cval );
			pos++;
		}
		map[iter] = r11;
		data[iter] = qgcomplex( fdata[pos], fdata[pos+1] ); pos+=2;
		iter++;
		assert( pos <= dsize ); 
	}  
	delete [] fdata;
	target.store_coefficients( map, data );
	return true;
}
// ---------------------------------------------------------------------------------------
template <int WIDTH>
bool read_static_var(const char * varname, map_container<REAL, WIDTH> & target, NcFile * ncfile){
	const int rwidth = 1; // require 2 elements to form a complex number
	int width, length = 0;
	if(! read_static_var_dims(WIDTH, width, length, rwidth, ncfile, varname) ) return false;
	// now read in variable
	NcVar * var = ncfile->get_var(varname);
	if( var == NULL || ! ( var->is_valid() ) ){
		fprintf(stderr, "Error: Variable '%s' not found in file.\n", varname);
		return false;
	}
	rank<INDEX, WIDTH> initializer;
	vector<rank<INDEX, WIDTH> > map( length, initializer );
	vector< REAL > data( length, 0 );
	int dsize = length*width;
	REAL * fdata = new REAL[dsize];
	if( ! var->get( fdata, length, width, 0, 0, 0 ) ){
		fprintf(stderr, "Error reading data from variable '%s'.\n", varname);
		delete [] fdata;
		return false;
	}
	// read data into index map etc ..
	int iter = 0, pos = 0;
	while(pos < dsize){
		rank<INDEX, WIDTH>  r11;
		for(int i = 0; i < WIDTH; i++){
			const REAL & cval = fdata[pos];
			if(cval < -127 || cval > 128){
				fprintf(stderr, "Error reading index map element '%i': Value exceeds INDEX limits.\n", iter);
				fprintf(stderr, "Error reading data from variable '%s'.\n", varname); delete [] fdata;
				return false;
			}
			r11[i] = static_cast<INDEX>( cval );
			pos++;
		}
		map[iter] = r11;
		data[iter] = fdata[pos]; pos++;
		iter++;
		assert( pos <= dsize ); 
	}  
	delete [] fdata;
	target.store_coefficients( map, data );
	return true;
}
// ---------------------------------------------------------------------------------------
bool read_static_data(nchandler & nc_container, const double & nchandle,
		map_container<REAL,3> &  Qcoeffs, map_container<REAL,3> & QIcoeffs){
	NcFile * ncfile = nc_container.getncfile(nchandle);
	if( ! read_static_var( "Q_operator", Qcoeffs, ncfile ) ){
		return false;
	}
	if( ! read_static_var( "QI_operator", QIcoeffs, ncfile ) ){
		return false;
	}
	return true;
};
// ---------------------------------------------------------------------------------------
bool read_static_data(nchandler & nc_container, const double & nchandle,
		static_container & Rs, static_container & Ss, map_container<REAL,6> & Acoeffs,
		map_container<REAL,3> &  Qcoeffs, map_container<REAL,3> & QIcoeffs){
	NcFile * ncfile = nc_container.getncfile(static_cast<int>(nchandle));
	if( ! read_static_var( "R_static_weights" , Rs, ncfile ) ){
		return false;
	}   
	if( ! read_static_var( "S_static_weights" , Ss, ncfile ) ){
		return false;
	} 
	if( ! read_static_var( "A_operator", Acoeffs, ncfile ) ){
		return false;
	}
	if( ! read_static_var( "Q_operator", Qcoeffs, ncfile ) ){
		return false;
	}
	if( ! read_static_var( "QI_operator", QIcoeffs, ncfile ) ){
		return false;
	}
	return true;
};
// ---------------------------------------------------------------------------------------
bool get_static_data(nchandler & nc_container, const double & nchandle, qgcomplex * & data,
		const char * name_real, const char * name_imag){
	NcFile * ncfile = nc_container.getncfile(static_cast<int>(nchandle));
	NcDim * sdim = ncfile->get_dim( "static_mapindex" ); // existence has already been checked for
	int map_size = sdim->size();
	if( map_size <= 0){
		fprintf(stderr, "Error: Static coefficient dimension has zero size.\n"); return false;
	}
	REAL * real;// = ncfile->get_var("R_real");
	REAL * imag;// = ncfile->get_var("R_imag");
	if( 
			! get_static_var( ncfile, real, map_size, name_real, sdim, ncREAL ) ||
			! get_static_var( ncfile, imag, map_size, name_imag, sdim, ncREAL )
	  ){
		fprintf(stderr, "Error reading static coefficients data.\n");
		return false;
	}
	data = new qgcomplex [ map_size ];
	for(int i = 0; i < map_size; i++ ){
		data[i] = qgcomplex ( static_cast<REAL>( real[i] ) , static_cast<REAL>( imag[i] ) );
	}
	delete [] real;
	delete [] imag;
	return true;
};
#endif
